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I The latest studies of massive star evolution indicate that an initially rapidly rotating 

, star with sufficiently low metallicity can produce a rapidly rotating, massive stellar core that 

could be a progenitor of long-soft gamma-ray bursts (LGRBs). Motivated by these studies, 
Ch ' we follow the collapse of a rapidly rotating massive stellar core to a 'collapsar' black hole (BH) 

^ I surrounded by a massive, hot accretion disk performing fully general relativistic simulations. 

' ' , We focus on the general relativistic dynamics of the collapse, and the relevant microphysics is 

00 ' treated in a qualitative manner. The simulations are performed until the system consisting of 

^N) ' the BH and the disk has relaxed to a quasi-stationary state. A novel mechanism found in this 

study is that strong shock waves are formed at the inner part of the disk after the formation 
of the BH. These shock waves propagate mainly along the rotation axis, heating the disk 
and sweeping materials around the rotational axis, and thereby forming a low density region. 
The temperature of the disk is high enough for copious neutrino emission. All these features 
^ , indicate that the direct formation of a rapidly rotating BH is a promising source of LGRBs 

^ ■ even in the absence of strong magnetic fields. 

a . 

I— I, §1. Introduction 

^ ■ The observed association of afterglows of long gamma-ray bursts (LGRBs) with 

'nJ" I star forming regions in galaxies P''^''^ such as the connections between GRB980425 

' ' and the Type Ib/c supernova SN1998bw,®''9 between GRB030329 and the Type Ic 

T-H ■ supernova SN2OO3dh,'Sl0® and between GRB060218 and the Type Ic supernova 

\^ ', SN2006aj|3'1^ (see also Ref . [TTI) for other suggested associations between LGRBs 

O I and core-collapse supernovae), suggested that at least some LGRBs are associated 

■ with supernovae and the core collapse of massive stars. Supported by these evidences, 

. . [ the collapsar scenario, in which the central engine of LGRBs is composed of a rotating 

' black hole (BH) surrounded by a massive, hot accretion disk,'I21'''I3''Hli ig currently 

^ ■ the most favored model for LGRBs. It requires the progenitors to be rotating rapidly 

^ I enough that an accretion disk can be formed around the BH. Also, the relativistic 

■ ■ ■ ' jets that are to produce the LGRB have to reach the stellar surface. This constrains 

the size of the progenitor; the progenitor should not have an extended hydrogen 
envelope. 'SI' '113 Thus, the progenitors of LGRBs are now believed to be rapidly 
rotating massive Wolf-Rayet (WR) stars. 

However, WR stars are known to be accompanied by strong stellar winds driven 
by radiation pressure which lead to a rapid spin-down of the cores of the WR stars.'^ 
Some author^^''!^ have proposed binary-evolution models for producing rapidly 
rotating progenitors of LGRBs. On the other hand, Yoon and Langer and Woosley 
and Heger have recently discoverecP^''^''^ that even a single star can fulfill the 
requirements of the collapsar scenario if it is initially rapidly rotating (> 50% of the 
Keplerian velocity at the equatorial surface) and of low metallicity {Z/Zq < 0.1). 
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The mass loss is suppressed by the low metallicity. ' The rapid rotation results in a 
short mixing timescale, which leads to a chemically homogeneous state throughout 
the hydrogen burning phaseP^ In this case, a single star can become a rapidly 
rotating WR star without losing the hydrogen envelope through the stellar wind, 
avoiding the red giant phase that otherwise would cause a significant decrease of the 
core angular momentum due to magnetic torques.'^''^ 

In addition to the above single-star evolutionary scenario to the LGRB pro- 
genitor, recent growing empirical evidence indicates that LGRBs may prefer a low 
metallicity environment. At least some LGRBs are likely to be formed 
from a rapidly rotating WR star which is born in a low metallicity environment and 
experiences chemically homogeneous evolution. 

Motivated by these latest studies, we perform axisymmetric simulations of rapidly 
rotating stellar core collapse to a BH in full general relativity. We explore the out- 
come of the collapse of very massive, rapidly rotating WR stars in the context of the 
collapsar scenario. As a first step toward more realistic simulations, microphysical 
processes are treated only in a qualitative manner, focusing on the general relativis- 
tic hydrodynamics of the collapse. Throughout this paper, c and G denote the speed 
of light and the gravitational constant, respectively. 

§2. Setting 

To model the core of WR stars just before the collapse, we adopt the polytropic 
equation of state (EOS) P = Kp^, with F = 4/3, taking account of degenerate 
electron and radiation pressures. Here is a constant. Rigid rotation is adopted as 
the rotation profile for simplicity. The maximum angular velocity (which is equal to 
the Kepler angular velocity of the equatorial surface) is assumed, since the progenitor 
(a WR star) of LGRBs is likely to be rapidly rotating. The central density and 
temperature of a pre-collapse iron core of massive metal-poor stars are pc ^ 10^ 
g/cm^ and Tc < 10^'' K (according to Umeda and Nomoto). We set the central 
density of the initial conditions to ~ 5 x 10^ g/cm^. Note that the chemical 
potential of the electrons at such a density is larger than the temperature T ~ lO^'' 
K, and hence the electrons are degenerate even at such a high temperature. 

According to the latest models of stellar evolution,!!^''^''^ rapid rotation results 
in a chemically homogeneous state and leads to a large GO core which could produce 
an iron core with mass > 2M0.'^ The CO core mass Mco can be approximately 10, 
17, 25Mq for models of initial mass 20, 30, AOMq, according to the result of Ref. l20p . 
Taking this fact into account, we choose models with a core mass of M ~ 3.5- 
4.5M0. This is achieved by setting K ^ 9-10.5 x 10^^ cm^/s^/g^/^. [Note that the 
mass is approximately given by 4.555(ir/G)3/2 for the T = 4/3 polytrope (see, e.g., 
Ref. [28]) ).] In this case, a BH is formed directly as a result of the collapse for any 
of the chosen EOSs (see the next paragraph). Because we have found qualitatively 
the same results for different mass models, we present the results for M = 4.2Mq 
in the following. In this model, the ratio of the rotational kinetic energy Trot to the 
gravitational potential energy W is Trot/|l^| ~ 0.0088. The nondimensional spin 
parameter q = cJ/GM'^ is ~ 0.98, where M and J denote the mass and angular 



Formation Mechanism of Collapsar Black Hole 



3 



Table I. The parameter set (A,-r2,Pnuc,/th), the maximum gravitational mass Mmax, and the time 
Iah at which the first apparent horizon is formed. 



Model 


A 




Pnuc(g/cm-^) 


Ah 


AAax (Mq) 


tAH (ms) 


A 


1.315 


2.6 


2.0x10'" 


1.315 


1.990 


230 


B 


1.32 


2.5 


1.9x10" 


1.32 


1.999 


276 


C 


1.325 


2.45 


2.0x10" 


1.325 


1.992 


358 



momentum of the core, respectively. 

For the evolution, we adopt a hybrid EOS (see, e.g., Ref. [29|) ) in which the 
pressure consists of the sum of cold and thermal parts, as P = Pcoid + -fthi where 



I Kip 1, p< pnuc, .X 
^2P , P > Pnuc, 



(Pnuc denoting the nuclear density) and Pth = (Ah — l)/o(e — ^coid)- Here, p is the rest- 
mass density, e and Ecoid are the total and cold parts of the specific internal energy, 
Fi, F2, and i~ih are constants, Ki = 5 x 10^^ in cgs units, and K2 = Kip^^~^^, 
respectively. By choosing a value of Fi slightly smaller than 4/3, the effect of the 
depletion of the degenerate electron pressure due to photo-dissociation and electron 
capture is qualitatively taken into account. The specific internal energy is set to the 
same value as in the case of the F = 4/3 polytropic EOS, i.e. e = 

The values of Fi are chosen as 1.315, 1.32, and 1.325 in this paper. For the 
case of ordinary supernova simulations (i.e. the collapse of a stellar core with a 
mass smaller than those in the present work), the value of Fi may be in the range 
1.30 < A < 1.33.123 For the collapse of cores with larger masses, it has been pointed 
out that the value of Fi tends to be larger than in the lower mass case, due to their 
larger entropies (see, e.g., Ref. [5CTj) ). For this reason, we believe that the values 
of Fi adopted in this paper may be appropriate. The values of F2 and Pnnc are 
chosen so that the maximum allowed gravitational mass of the cold neutron star 
is Mjnax ~ 2.OM0, which is approximately equal to the highest pulsar mass ever 
measured,!^ i.e., M = 2.1 ± O.2M0. We set Fth = A for simphcity. With such 
small values, the low efficiency of the shock heating due to the energy loss through 
photo-dissociation and neutrino emission is qualitatively taking in to account . The 
adopted values for the parameter set (AjAi/fnucjAh) are listed in Table I. 

The temperature T is estimated from e by solving the equation 

e = egas + erad + £iy + Ecold , (2-2) 

where Egas, eradj and denote the specific internal energy of the gas, the radiation, 
and the neutrinos and are functions of p and T. We use the same forms for these 
quantities as in Ref. |32]) , setting 

_ 3pkT 1 + 3Xnuc „x 
2mp 4 ' ^^'"^^ 

erad = ^aF\ (2-4) 
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= ^ar40(r-r„it), (2-5) 

o 

he /3^y/3 

where k is the Boltzmann constant, h Planck's constant, rup the nucleon mass, Xnuc 
the mass fraction of free nucleons, a the radiation constant, and the electron 
fraction. Further, denotes the Heaviside step function; i.e., for an optical depth r 
larger than the critical optical depth Tcrit ~ 2/3, we assume that the neutrinos are 
opaque with respect to nucleons and electrons and take into account their energy 
densities. (The definition of r is given in §3.) 

Simulations were performed using an axisymmetric, full general relativistic codeW^ 
The so-called simple excision techniqud^S was implemented to follow the evolution 
of the BH spacetime in a stable manner. (For recent calculations, see, e.g., Refs.EJj) 
and i36i) .) The regridding techniqud^ was adopted during the collapse. The grid size 
and grid spacing in the final regridding were (3300, 3300) for (cj, z) and Ax w Af/20, 
where lo = y^x^ + y-^ is the radius in cylindrical coordinates. 

§3. Results 

In the early stage, during which the central density of collapsing core is smaller 
than the nuclear density, the collapse proceeds in an approximately homologous man- 
ner. Then, as the central density approaches the nuclear density, the collapse around 
the equatorial plane is decelerated, due to the centrifugal force, while the collapse 
along the rotational axis is relatively accelerated to form a flattened structure. Since 
the mass of the stellar core is much larger than the maximum allowed mass for the 
given EOS, a BH is formed directly. The formation of the BH is ascertained by find- 
ing the apparent horizon. After the apparent horizon is formed (for the formation 
time, see Table I), a centrifugally supported, thin accretion disk is formed around 
the BH (see the top left panel of Fig. 1): The disk rotates at approximately the Ke- 
plerian velocity, (see the bottom panel of Fig. 2). The evolution of the irreducible 
mass of the BH, Mirr = ^/A/IQ-k, and the disk mass Mdisk are shown in Fig. 3. Here 
A is the area of the apparent horizon. 

After the disk formation, shocks are formed at the inner part of the disk, 
converting the kinetic energy of the infall into thermal energy (See Ref. [38]) for 
discussion of a similar phenomenon.) The gravitational potential energy released 
is ^ ~ GMBHMdisk/^isco ~ 4-5x10^2 g^gs, where Mbh, Mdisk ~ O.IMq, and 
RiSCO ~ 4-5Gc~^Mbh are the black hole mass, the disk mass, and the radius of the 
innermost stable circular orbit around the BH, respectively. As the thermal energy 
is stored, the disk height H increases. For H < -RiscO) H in the vicinity of the BH 
is approximately determined by the balance relation 

-Pdisk — -Pram GMbrPsH 



-"^isco 



(3-1) 



where the left-hand and right-hand sides represent the pressure gradient and the 
gravity of the BH. The quantities ps, -Pdiskj and Pram are the disk density near the 
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Fig. 1. Snapshots of density contours and velocity fields of the central region of the core in the x-z 
plane at t ~ 280.82ms (the top left panel), 281.60ms (the top right panel), 282.25ms (the middle 
left panel), and 290.14ms (the middle right panel). Snapshots of density contours (the bottom 
left panel) and temperature contours (the bottom right panel) in the x-z plane at 290.14 ms. 
All figures are for model B. 



surface, the pressure inside the disk, and the ram pressure of the infalUng matter, 
respectively. Equation (|3-ip provides the order of magnitude of the pressure as 
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Fig. 2. Various quantities characterizing the disk at 290.14 ms for model B. The top panel plots 
the surface density Sn = E /IQ^"^ g/ca? (solid curve) and the optical depth r s» k^S (dotted 
curve). The middle panel plots the temperature evaluated in the equatorial plane. The bottom 
panel plots the angular velocity Q (solid line) and the Keplerian angular velocity (dotted curve) . 
The vertical dashed lines at J? w 19 km indicate the expected location of the marginally stable 
circular orbit around the BH. 



The ram pressure can be expressed as 

Pram ^ Pivj ~ lO^o f 3 ) dyn/cm^ (3-3) 

VlO^^g/cm / 

where pf and Vf ~ (2GMbh/-Risco)^^^ ~ 0.4-0. 5c are the density and velocity of the 
inf ailing matter. 

The density (paisk) and the temperature (Tdisk) inside the disk eventually increase 
to ~ 10^^ g/cvc? and ~ lO^-*^ K (and hence, Paisk ~ lO^"*^ dyn/cm^), while the ram 
pressure (i-*ram) decreases to < ClPjisk {"^ -Pdisk); since the density of the infalling 
matter (pf) decreases to < lO^'^ g/cm^. Consequently, H increases to ~ Pisco foi' 
Ps ~ 10^^ g/cw?. Then, the approximate balance relation becomes 

(Pdisk - Pram) ~ ■ (3-4) 

Since the binding due to the gravitational force decreases as H increases, the disk 
expands preponderantly in the z-direction, forming strong shock waves (see the top 
right and middle left panels of Fig. 1). The propagation speed of the shock waves 
is K, 0.5c, i.e., mildly relativistic. As the shock waves propagate, the mass accretion 
rate to the BH significantly decreases (see Fig. 3). After the shock propagation, a 
low density funnel region is formed around the rotational axis (see the middle right 
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Fig. 3. Evolutions of the irreducible mass Ali^ and the disk mass i\fdisk- The long dashed, solid, 
and dashed lines represent the results for models A, B, and C, respectively. The horizontal line 
at M 4.2M0 indicates the total mass. 



and bottom left panels of Fig. 1). In addition, the funnel region is surrounded by a 
dense, hot wall (see the bottom right panel of Fig. 1). The temperature in the funnel 
increases to approximately the same level as that in the wall through shock heating, 
and hence such a low density funnel can be supported by the thermal pressure. 

Figure 2 displays the surface density, the optical depth, the temperature, and 
the angular velocity of the disk at t = 290.14 ms. Due to the shock heating, the 
temperature increases to T ~ 4-9 x 10^*^ K, which is about a factor of 4-5 larger than 
that in the unshocked regions (see the middle panel of Fig. 2). The surface density 
around the rotational axis is S « 2-6 xlO^^g/cm^, while it is T ~ 10^^ g/cm^ in the 
disk. This reflects the formation of the funnel structure formation (see the top panel 
of Fig. 2). Assuming that neutrino opacity is = 7 x 10~^''(T/10^^ K)^ (see, e.g., 
Di Matteo et al. 2002), we compute the optical depth r = k^^- Then we find that 
the disk is optically thick for a radius of f« 20-130 km (see the top panel of Fig. 
2). The subsequent evolution of the hot disk will be determined by the duration of 
the infall and the viscous timescale, as previously discussed (see, e.g., Ref. [39|) ). 

Successful LGRB formation requires at least the following three ingredients: 1. 
formation of an energetic fireball which yields highly relativistic outflows; 2. colli- 
mation of the relativistic outflows; 3. successful penetration of the outflows through 
a surrounding stellar mantle. As discussed in the following, all these conditions are 
likely to be satisfied in our model, even in the absence of a strong magnetic field, 
which is often assumed (see, e.g., Refs. HOj) and HI])). 
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In our models, a highly relativistic fireball will be produced by the neutrino pair 
{vv) annihilation around the rotational axis. Here, let us approximately estimate 
the energy deposition rate Eyp through uv annihilation. In the diffusion limit, the 
neutrino ffux is given by 

Fu ~ — ^, (3-5) 

where Ni, = 3 is the number of neutrino species and a is the Stefan-Boltzmann 
constant. The neutrino luminosity from an optically thick disk is then 

« 27ri?Lki^. ~ 5 X lO^^ ^-i {^^^^ ergs/s (3-6) 

where Tn = T/IQ^^ K and Un = 17/10^ g/cm^. According to the results of 
Setiawan et al.,H3 the expected energy deposition rate through the z^z/ annihilation 
would then be E^p ~ 10^^ ergs/s. The low density funnel region above the BH will 
be a favorable place for an efficient production of ee"*" pairs through z/P annihilation. 

The formation of a dense, hot wall surrounding the funnel will play an impor- 
tant role in collimating relativistic outflows. In the absence of this wall, energetic 
outflows fail to be collimated, as the outflows travel through the stellar mantle. The 
shock waves formed at the birth of the BH are likely to sweep the matter along 
the rotational axis, reaching the stellar surface before the main relativistic outflows, 
which will be produced in the accretion phase (e.g., Ref. [39])), are driven. 

For larger values of l-Ti— 4/3|, the time at the onset of the shock wave propagation 
is delayed, because the rate of increase of the disk pressure as well as the rate of 
decrease of the ram pressure are smaller, reflecting the less homologous nature of the 
collapse. This also results in the fact that it takes a longer time for the BH to relax 
to a stationary state (see Fig. 3). However this 'time delay' is at most ~ 10 ms for 
Fi > 1.315, and the features of the funnel structure formation are qualitatively the 
same for any of the adopted EOSs. 

The thermal energy is stored in the inner part of the accretion disk if the heating 
timescale is shorter than the timescale of the neutrino cooling. The condition is 
a{GMm/R) ^ Ly (~ 10^^ ergs/s), where a < 1 is the effective conversion efficiency 
of the kinetic energy to the thermal energy through the shocks and m is the mass 
accretion rate. This gives m S> LyR/ aGM ~ a~^A'lQ s"^. According to our results, 
m > IOM0 s"^ at least for {t - iAn) < 10 ms. Thus, the neutrino cooling does 
not play an important role in the thermal-energy-store phase, unless the conversion 
efficiency is extremely low, i.e., a <C 0.1. To summarize, our results are universal, at 
least qualitatively, for any of the EOSs and microphysical processes. 

We also performed simulations for 0.8 ^ q 1.2, as well as for weakly differen- 
tially rotating cases and found essentially the same results for all these parameter 
values. For larger values of q ~ 1.2, more collimated shock waves are formed. For 
smaller values of g < 0.8, a massive disk, which is essential for shock formation, is 
not formed. For sufficiently large values of g > 1.2, on the other hand, a BH is not 
promptly formed, due to the effect of the centrifugal force. In this BH may 

be formed after a sufficient amount of angular momentum is transported. 
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The black hole excision technique enables us to continue the simulation until 
the BH has relaxed to a quasi-stationary state. As Fig. 3 shows, a seed BH with 
mass Mirr ~ 2AMq is born at first for any of the EOSs. This value is approximately 
determined by the maximum allowed mass of the rigidly rotating configuration for 
the adopted EOSs. The figure also reveals that approximately 95% of the total 
mass is eventually swallowed by the BH for any of the adopted EOSs. Since the 
total angular momentum is conserved, the angular momentum of the BH can be 
indirectly estimated from that of the disk, which is ~ 20% of the total angular 
momentum. Thus, the final spin parameter of the BH is estimated as « 0.8. 

§4. Summary 

Motivated by recently developed single-star evolutionary models of LGRB pro- 
genitors, we performed fully general relativistic simulations for rapidly rotating stel- 
lar core collapse to a BH. We found that a BH is directly formed as a result of the 
collapse of a sufficiently massive progenitor. Soon after the BH formation, a disk of 
density ^ 10^^ g/cm^ is formed around the BH. The subsequent infall of matter from 
outside produces shocks at the surface of the disk, and thermal energy is generated, 
heating the disk. The thermal pressure eventually reaches ~ 10^^ dyn/cm^, which 
is much larger than the ram pressure of the infalling matter, and hence it is used to 
sustain the vertical gravitational force. Because of the increasing thermal energy, the 
disk and shock front gradually expand, and when the scale height is comparable to 
the disk radius, the shock front starts to expand with mildly relativistic speed. The 
strong shock waves sweep matter around the rotational axis and heat up the disk 
matter to ~ 10^^ K. Then, a low-density {p < 10^ g/cm^) funnel region surrounded 
by a hot, dense wall is formed. Formation of such a structure in advance of the 
subsequent quasi-stationary accretion from the hot disk and the resulting relativistic 
outflows will be quite favorable for producing the fireball and LGRBs. 

In the present simulation, we did not incorporate magnetic fields. It should 
be noted that even in the absence of magnetic stress, funnel structure, which will 
be essential for producing collimated jets, is formed. However, magnetic fields can 
further promote the formation of such a funnel if the field strength is larger than 
~ 10^^ G (i.e., the magnetic pressure is larger than ~ 10^^ dyn/cm^, which is 
comparable to the ram pressure of the infalling matter). Study of magnetic field 
effects is left as a future project. 

Although our treatments of the EOS and microphysical processes in the present 
work are approximate, we believe that the qualitative features of the BH formation 
and subsequent shock formation found here will be universal. To confirm this process 
more rigorously, more detailed quantitative studies with a realistic EOS, including 
relevant microphysics, is needed. Such a study is in progress, and the results will be 
reported in the near future. 
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